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Abstract 

The equation of state and composition of the inner crust of neutron stars at zero temperature 
are calculated, using the T = version of the TETFSI (temperature-dependent extended Thomas- 
Fermi plus Strutinsky integral) method, for each of a family of three functionals based on Skyrme- 
type forces BSkl9, BSk20 and BSk21, which are characterized by different degrees of symmetry- 
energy stiffness, and also for the SLy4 functional. We also solve the Tolman-Oppenheimer-Volkoff 
equations to calculate the distribution of mass within the inner crust. Qualitatively similar results 
are found for all four functionals, and in particular the number of protons per Wigner-Seitz cell 
is in all cases equal to 40 throughout the inner crust. 

PACS numbers: 04.40.Dg, 21.10.Dr, 21.30.-x, 21.60. Jz, 26.60.Gj, 26.60.Kp 
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I. INTRODUCTION 



We recall that three distinct regions can be recognized in a neutron star: a locally homoge- 
neous core and two concentric shells characterized by different inhomogeneous phases [l[Q]- 
The outermost of the shells, the "outer crust" , consists of an electrically neutral lattice of 
nuclei and electrons. At the surface of the star only nuclei that are stable under natural 
terrestrial conditions are found (in fact, under the assumption of "cold catalyzed matter", 
i.e., nuclear and beta equilibrium at temperature T = 0, only 56 Fe will be found), but on 
moving towards the interior the increasing density leads to the appearance of nuclei that are 
more and more neutron rich, until at a mean local density n of around 2.5 xlO -4 nucleons 
fm -3 (4.2 x 10 11 g cm -3 ) neutron drip sets in. This marks the transition to the "inner crust" , 
an inhomogeneous assembly of neutron-proton clusters and unbound neutrons, neutralized 
by an essentially uniform electron gas. By the point where the mean density has risen to 
about two thirds of the density no of symmetric (homogeneous) nuclear matter (SNM) at 
equilibrium, the inhomogeneities have been smoothed out and we enter the core of the star. 
The homogeneous medium of which the core is comprised is known as "neutron-star matter" 
(N*M), and is made up primarily of neutrons, with a small admixture of protons neutral- 
ized by electrons (and muons at densities above n ~ 0.12 fm -3 ). Closer to the center, other 
particles such as hyperons might appear. 

In this paper we continue our calculations of the different regions of neutron stars with a 
family of three Skyrme-type functionals, BSkl9, BSk20 and BSk21, that we have constructed 
specifically to provide a unified approach not only to the structure of the different regions of 
neutron stars but also to other phenomena associated with the birth and death of neutron 
stars, e.g., supernova-core collapse, the r-process of nucleosynthesis in the neutrino-driven 
wind, and nucleosynthesis via the decompression of neutron-star matter j^|. These three 
functionals are all based on effectives forces with the generalized Skyrme form 

Vij = * (1 + xoP^irij) + + xiP*)jz [p 2 l3 S(rij) + 5(r^)pJ] 

+ t 2 (l + x 2 P a )-^p ij .8{r ij )p ij + -t 3 (l + x 3 P a ) n(r) a 5(r i:j ) 

+ -U(l +x 4 P CT )^ [pln{rf 5{ riJ ) + 5{r l0 )n{rf pi] 

+ t 5 (l + x 5 P a )-2P ij .n{r) 1 5(r lj )p ij + Xw (tTi + aj) ■ p id x S^p^ , (1) 
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where = Ti —Tj, r = (r, +rj)/2, = —ih(Vi — Vj)/2 (this is the relative momentum), 
P CT is the two-body spin-exchange operator, and n(r) = n n (r) + n p (r) is the total local 
density, n n (r) and n p (r) being the neutron and proton densities, respectively. The £ 4 and 
t 5 terms here are unconventional, being density-dependent generalizations of the t\ and £2 
terms, respectively. 

The parameters of this form of force were determined primarily by fitting measured nu- 
clear masses, which were calculated with the Hartree-Fock-Bogoliubov (HFB) method. For 
this it was necessary to supplement the Skyrme forces with a microscopic contact pairing 
force, phenomenological Wigner terms and correction terms for the spurious collective en- 
ergy. However, in fitting the mass data we simultaneously constrained the Skyrme force to 
fit the zero-temperature equation of state (EOS) of homogeneous neutron matter (NeuM), 
as determined by many-body calculations with realistic two- and three-nucleon forces; the 
strength of the pairing force at each point in the nucleus in question was likewise calculated 
analytically so as to reproduce the x Sq pairing gaps of homogeneous nuclear matter of the 
appropriate density and charge asymmetry |4j. Actually, several realistic calculations of 
the EOS of NeuM have been made, and while they all agree very closely at nuclear and 
subnuclear densities, at the much higher densities that can be encountered towards the cen- 
ter of neutron stars they differ greatly in the stiffness, i.e., the density dependence, of the 
symmetry energy that they predict, and there are very few data, either observational or 
experimental, to discriminate between the different possibilities. It is in this way that we 
arrived at the three different functionals of this paper: BSkl9 corresponds to the softest 
EOS of NeuM known to us, BSk21 to the stiffest, while BSk20 has intermediate symmetry 
stiffness, as seen in Fig. 1 of Ref. {3]. On the other hand, Fig. [I] of the present paper shows 
that in NeuM the three functionals are very close to each other at the subnuclear densities 
relevant to neutron-star crusts. For a further discussion of this point see Ref. {3J, where it 
will be seen in particular that a value of 30 MeV was imposed on the symmetry coefficient 
J for all three functionals. It will also be seen there that the values of the density-symmetry 
coefficient L, which measures the stiffness of the symmetry energy at the equilbrium density 
Uq, are all very similar. 

Furthermore, we imposed on these functionals the supplementary constraints of i) elim- 
inating all unphysical instabilities in nuclear matter for all densities up to the maximum 
found in neutron stars (these functionals are also stable at the finite temperatures encoun- 
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FIG. 1: (Color online.) Neutron-matter EOSs (internal energy per nucleon e as a function of 
density n) for forces BSkl9 - 21 and SLy4 at subnuclear densities and zero temperature. 

tered in supernova cores s|) ii) obtaining a qualitatively realistic distribution of the potential 
energy among the four spin-isospin channels in nuclear matter hi) ensuring that the isovector 
effective mass is smaller than the isoscalar effective mass, as indicated by both experiment 
and many-body calculations. 

The introduction of the unconventional terms in £4 and t§ allowed us to satisfy all these 
constraints and at the same time fit the 2149 measured masses of nuclei with N and Z > 
8 given in the 2003 AME (Atomic Mass Evaluation) [6] with an rms deviation as low as 
0.58 MeV for all three models, i.e., for all three options for the high-density variation of the 
symmetry energy. For all three of these functionals complete mass tables (labeled HFB-19, 
HFB-20 and HFB-21, respectively) were constructed, going from one drip line to the other. 
The reliability of the predictions that these models make for experimentally inaccessible 
neutron-rich nuclei is all the greater for the constraints to neutron matter imposed on their 
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underlying forces, and it wa 8 thus particularly appropriate to use these mass models in our 

earlier study of the outer crust of neutron stars |7|]. As for the homogeneous core, the T 

= EOSs of N*M for our forces BSkl9, BSk20 and BSk21 have already been published in 

I I 

the original paper presenting these forces [3J. This leaves just the inner crust to be dealt 
with, and our main concern in this paper is to calculate for this region the EOS and the 
composition as a function of density with each of our three functionals. 

We shall also perform inner-crust calculations with the functional SLy4 jg], since like 
our own functionals it is designed for finite-nucleus HFB calculations, and is intended for 
neutron-star studies, being subject to a neutron-matter constraint. However, it has the 
conventional Skyrme form and thus, having fewer parameters, is far less flexible than our 
own functionals. Thus SLy4 was fitted to only six nuclear masses; moreover, three of these 
nuclei had N = Z (even), and since no Wigner term was included in the model the symmetry 
energy must inevitably be too large. (In particular, the symmetry coefficient J for this 
functional is 32 MeV, while we have found that the optimal value for the conventional 
form of Skyrme functional when all the mass data are fitted without any neutron-matter or 



other constraint is 28 MeV 



9].) The excessive symmetry energy might explain why the rms 



lOt ] : note that only even-even nuclei 



deviation from the mass data is quite large, 5.1 MeV 
were considered in that calculation. 

Given that all four functionals were fitted to masses with the HFB method, it might 
seem appropriate to use this method for the inner-crust calculations as well. Now the latter 
calculations have been generally performed within the framework of the spherical Wigner- 
Seitz (WS) approximation, as in the pioneer HF calculations of Negele and Vautherin 111 ], 
in order to avoid computer times grossly in excess of those for isolated-nucleus calculations. 
But an inevitable consequence of the WS approximation is to introduce shell effects in the 
spectrum of unbound neutron states, which dominate the properties of the inner crust. Such 
shell effects are to a large extent spurious, since in reality the unbound neutron states form 



a quasi- continuum. This difficulty is analyzed in detail in Refs. 12 



the latter reference 



showing that the error thereby introduced in the energy per nucleon cannot easily be reduced 
below 50 keV, which is incompatible with a reliable calculation of the composition of the 
inner crust; for a very recent discussion of the problem see Gril 
years 3D calculations have been carried out by several groups 



15 



et al. 1 141. In the last few 



-Il7j]. However, not only 



does this sort of calculation require computer times that are quite impractical for extensive 
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astrophysical calculations but the use of a cubic box with periodic boundary conditions can 
still lead to spurious neutron shell effects (see, for example, Section C.2 in Ref. 

In view of these problems it is not surprising that a more popular approach to the 
calculation of the inner crust has been to use the much simpler compressible liquid-drop 
model (CLDM); a typical such calculation is that of Ref. 18) . Within each Wigner-Seitz 
cell this method makes a clear separation of nuclear matter into two distinct homogeneous 
phases, the densities of which are free parameters of the model. The bulk properties of the 
two phases are calculated microscopically using the adopted functional, as are the surface 
properties (preferably including curvature corrections) of the interface between them. A 
more realistic treatment of spatial inhomogeneities is to employ semi-classical methods such 
as the Thomas- Fermi (TF) approximation, as for instance in Ref. 19j. However, both 
the CLDM and TF methods are otherwise purely macroscopic, and in particular have no 
quantum shell corrections at all. 

The so-called TETFSI method (temperature-dependent extended Thomas-Fermi plus 
Strutinsky integral) of Onsi et al. [20j, which we adopt here, is a computationally very 
fast approximation to the full finite-temperature HF method. This method was originally 
developed for calculating the EOS of the dense matter found in supernova cores [2l|. But in 
this work we will be using just the zero-temperature limit. The TETFSI method, like the 
TF method, allows for a continuous variation of the density of nuclear matter within each 
WS cell, without any artificial separation into two distinct phases. However it is expected 
to provide a much better description of nuclear clusters than the TF method because the 
semi-classical expressions for the kinetic-energy and spin current densities include density- 
gradient terms up to the fourth order. Most importantly, proton shell corrections are added 
perturbatively, but we avoid the difficulty of spuriously large values for the neutron shell 



corrections noted above by not calculating them at al 



much smaller than the proton shell corrections 13 



22 



; in any case they are known to be 



Our method is described in detail in Ref. 20], but we summarize it here in Section II. The 
results for the zero-temperature composition and EOS of the inner crust are presented in 
Section III, along with an examination of the extent to which continuity holds at the interface 
with the outer crust. This section also discusses the transition between the inner crust and 
the liquid core. In Section IV we examine the solutions to the Tolman-Oppenheimer-Volkoff 
(TOV) equations 24], |25) in order to determine the distribution of mass within the inner 
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crust. Our conclusions are summarized in Section V. 



II. TETFSI MODEL OF INNER CRUST 



To summarize the main features of the TETFSI method 20j, we note first that it models 
the inhomogeneous medium by spherical WS cells, with the spherically symmetric neutron 
and proton density distributions being parametrized according to 

n q {r) = n Bq + n Aq f q (r) , (2) 

in which, with q = n or p, n Bq is a constant background term, while 

h(r) = r~, A — s ; r > ( 3 ) 



Cq —R\ -i I f r—C, 



l + exp ^ -1 exp 



R 

In this "damped" form of the usual simple Fermi profile all density derivatives vanish at 
the surface of the cell, thereby ensuring a smooth matching of the nucleonic distributions 
between adjacent cells, and satisfying certain necessary conditions discussed below. It is 
particularly to be noted that with this parametrization of the density there is no arbitrary 
separation into liquid and gaseous phases within the WS cell. However, if this were what 
is energetically favored in reality, it would automatically be taken into account through the 
small values of the diffusenesses a q that would emerge. 

In order to determine the composition and the EOS of the inner crust one should in 
principle minimize at constant pressure the Gibbs free energy g per nucleon with respect 
to all the parameters of the WS cell. This is the procedure that we adopted in Ref. [7J for 
the outer crust, but for the inner crust the computation would be extremely heavy. Instead, 
we minimize rather the total Helmholtz free energy / per nucleon at constant mean density 
n with respect to the same parameters, showing in Appendix |A] that the error thereby 
introduced is quite negligible. Since the present work is limited to T = it is the internal 
energy per nucleon e that is minimized (/ = e — Ts, where s is the entropy per nucleon). 

To enumerate the minimizing parameters of the WS cell, we note first that the cell 
radius R will be determined, for the given n, by the total number A of nucleons in the cell. 
Then with the number of protons Z and the number of neutrons N in the cell specified 
(Z + N = A), only three of the four remaining cell parameters appearing in Eqs. (j2]) and (J3]) 
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for each charge-type of nucleoli will be independent. Thus, including Z and N, there will 
be eight parameters with respect to which the energy e must be minimized. Identifying the 
different contributions to e, we write 



+ e e + e c -Y e Q n> p , (4) 



and now discuss briefly each term. 
The nuclear term is 



e nuc = ^ / r^ F (r)dr + e^ , (5) 

n J cell 

where £fj,T F ( r ) is the ETF approximation to the energy density £sk y (?") given by Eq. (A3) 
of Ref. 26J for the generalized Skyrme force (JTJ (the formalism of Ref. 20] is limited to 
conventional Skyrme forces). All terms in Ssky are functions of the number densities n q (r), 
the kinetic-energy densities r q (r) and the spin-current densities J q {r). The ETF method 
approximates these last two densities as functions of the number densities n q (r) and their 
first four derivatives. However, as far as the total ETF energy is concerned, it is shown in 
App. A of Brack et al. [27] that the third- and fourth-order derivatives of the density can be 
eliminated by partial integration over the region of interaction, provided certain boundary 
conditions are satisfied on the bounding surface of this region. In the finite-nucleus case of 
Ref. 27J] the bounding surface can be taken to lie at infinity, in which case the necessary 
boundary conditions are trivially easy to satisfy. In the present case the bounding surface 
is the surface of the WS cell, and, given the fact that in general the density does not vanish 
on this surface, the necessary conditions are that the first three derivatives of the density 
must vanish there. These conditions are satisfied automatically for the distribution ([3]), and 
it is in this "integrated" form that we have implemented the ETF method. Note that for all 
four functionals, i.e., BSkl9-21 and SLy4, we omit the spin-current terms in J 2 , since this 
is the way these functionals were fitted (see Ref. |f| for a discussion of the implications of 
these terms for spin and isospin stability). 

With the ETF approximations for r q (r) and J q (r) being semi-classical all shell effects in 
£sky(f) are lost. The second term on the right-hand side of Eq. ((5]) represents our attempt to 
restore the proton shell corrections perturbatively using the SI (Strutinsky integral) method, 



as described in Ref. 



201 ]. As explained in Section I, we do not calculate neutron shell 



corrections in the inner crust; for a fuller discussion of this point see Section I of Ref. [20], 



where we conclude that because of the problems with neutrons, the ETFSI method is better 
adapted to a WS approach than is the HFB (or HF-BCS) method. On the other hand, we 
do not include pairing at the present stage of our calculations. This should have very little 
impact on the EOS, but it might have implications for the composition. 

The term e e on the right-hand side of Eq. (T4]) denotes the kinetic energy per nucleon of the 
electrons. In dense, cold, neutron star crust, electron-ch arg e screening effects are negligible 



and the electron density n e = n p is essentially uniform 



28 



29j . The energy e e can thus be 



calcu 
Ref. 



ated straightforwardly by expressions given in Section 24 of Cox and Giuli [30|, as in 



20|. 



The third term on the right-hand side of Eq. (J4J) denotes the total Coulomb energy per 
nucleon. It is calculated according to Eq. (3.4) of Ref. [20], except that there are the following 
changes to the exchange part, a) The proton exchange energy is set equal to zero for the 
three BSk functionals; this is a device that we have successfully adopted in all our recent 
models, beginning with BSkl5 [31]. and it can be interpreted as compensating for neglected 
effects such as Coulomb correlations, charge-symmetry breaking of the nuclear forces, and 
vacuum polarization, b) The electron exchange energy, which has the nonrelativistic form 
in Eq. (3.4) of Ref. |20| . is multiplied by a factor of -1/2, as appropriate for extremely 
relativistic particles [32]. 

The last term on the right-hand side of Eq. (J3J), in which Q n ^ is the beta-decay energy of 
the neutron (0.782 MeV) and Y e = Z/A, takes account of the neutron-proton mass difference 
(we drop a constant term M n c 2 ). 

Minimization of e with respect to the eight available parameters is performed by means of 
the CERN routine MINUIT. Actually, we found it necessary to exclude the shell correction 
term e p sh from this minimization, and then to add it later to what is really just the optimal 
ETF part of the energy. Otherwise, the minimization routine will tend to seek large negative 
values of the shell corrections, in violation of the essentially perturbative character of the SI 
method. In practice, we performed the minimization for different fixed values of Z, thereby 
reducing the number of free variational parameters to seven. However, even with this reduced 
number of parameters MINUIT occasionally failed to find a correctly converged minimum. 
This problem could often, but not always, be avoided by adjusting the initial values for the 
parameters. When this procedure failed solutions could always be found, provided we are 
not too close to the interface with the core, by a slight shift in the value of n; for this reason 
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our grid of values of n is irregular. However, above a certain value of n we were unable 
to find any solutions at all when MINUIT minimizes with respect to seven variables. We 
attribute the failure of our code to find well defined minima before true homogeneity has 
been reached to the energy minimum being very flat, with the result that MINUIT is unable 
to pick out one configuration among a very wide range of possibilities. We found, however, 
that we could still find well defined minima in this region by reducing the number of free 
variables in MINUIT to three, ri\ n ,n\ p and N, and minimizing for a large number of fixed 
values of the other five parameters; clearly, for a given level of accuracy this procedure will 
require much more computation time than when MINUIT minimizes on seven variables. 

It should be noted that at all densities the number of neutrons N in the WS cell is taken 
as one of the minimizing variables in MINUIT and hence is treated as a continuous variable, 
rather than being discretized to integral values. Even though the total number of neutrons 
in the crustal layer is, of course, integral, the notion of a fractional number of neutrons per 
WS cell corresponds, in fact, to the physical reality, since the neutrons are delocalized. 

Normally we would expect positive values of the constants n\ n and n\ p to emerge from 
the minimization, the cluster then representing a "droplet". However, there have been 
indications that towards the interface with the core the clusters may take several 

other forms. Most of these "pasta" configurations, such as slabs, tubes and rods, cannot be 
handled by our code, which is restricted to spherical shapes, but another of these possibilities, 
spherical bubbles, could in principle emerge from the minimization with our code, since they 
correspond simply to negative values of n\ q . We return to this possibility in Section Ull CI 

The pressure P corresponding to any given value of n is calculated by evaluating a simple 
analytic expression, as described in Appendix [B] This is more reliable and computationally 
much faster than the numerical differentiation of e used in Ref. 1 201 . 




III. COMPOSITION AND EQUATION OF STATE OF INNER CRUST 
A. Generalities 

Following the methods described in the previous section, for each of our three functionals 
and SLy4 we minimized the internal energy per nucleon e at temperature T = for more 
than a hundred different densities n between the drip point and 0.1 fm -3 . At this upper 
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limit our density distributions have become effectively homogeneous, as will be discussed in 
more detail in Section IIII CI 

For all values of n up to 0.06 fm -3 the optimal value of the number of protons Z per 
Wigner-Seitz cell was found to be 40, for all four functionals. However, at higher densities, 
as homogeneity is approached, the minimized energy becomes increasingly insensitive to Z. 
The preference for Z = 40 in the case of these four functionals is somewhat fortuitous, given 
that for some of our older functionals different values were found. For example, with the 
functional BSkl4 used in Ref. 



20] it was found that Z could take any of the values, 20, 40 
and 50, according to the density. It is remarkable that these familiar finite-nucleus magic 
proton numbers should persist in the highly neutron-rich environment beyond the drip line, 
especially in view of the presence of electrons, which will have the effect of significantly 
reducing Coulomb effects. 

For the specific case of functional BSkl9, reference to Fig. [2] shows both the role of shell 
effects and the overall trends imposed by the ETF part of the calculation. For both of the 
extreme densities shown here the ETF minimum lies close to Z = 40, and the shell effects 
simply reinforce this preference. However, the energy difference per nucleon A e between Z 
= 40 and Z =50 is very small: about 10 keV at the drip density and 5 keV at n — 0.06 
fm -3 (note the different energy scales of the two panels). It is easy to see how with even an 
only very slightly different functional a quite different T = composition could be found, 
as a result of changes in either the shell effects or the macroscopic ETF part of the energy 
(or both). 

Since the functionals BSkl9, BSk20 and BSk21 give better and wider data fits than all 
our earlier functionals, and have a better theoretical base as well, we believe our prediction 
of Z = 40 at all densities in the inner crust to be more credible than our earlier predictions, 
but the need for caution is evident. For example, taking pairing into account might well 
shift the favored value of Z away from 40. In any case, in a real neutron star a fairly wide 
range of values of Z can be expected at any point in the inner crust because of the finite 
temperature. 

The optimal values of A are plotted as a function of the density in Fig. [3J similarly, 
Figs. H] and [5] show the variation of e and the pressure P, respectively; these two figures 
show the densities n^ni °f transition between the inner crust and the core, as calculated in 



Section HHn 
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FIG. 2: Variation of ETFSI energy e per nucleon as a function of Z for functional BSkl9 with 
N optimized for each value of Z; dotted curve represents ETF approximation. Upper panel: 
n = 2.63 x 10 -4 nucleons fm -3 (drip density); lower panel: n = 0.06 nucleons fm -3 . 

No essential differences will be perceived between any of these four functionals, as far as 
the inner crust is concerned, although BSk21 is seen in Fig. [5] to have a somewhat softer 
EOS (in contrast to a much stiffer EOS at high density). These features can be related 
to the behavior of the respective functionals in homogeneous NeuM at inner-crust densities 
(see Fig. [p. Since SLy4 gives a much worse mass fit than do any of the BSk functionals, 
one might have expected that it would represent less well the presence of inhomogeneities 
and of protons, and thus give significantly different results in the inner crust, but this turns 
out not to be the case. Furthermore, the higher value of the symmetry coefficient J in the 
case of SLy4 (32 MeV, as opposed to 30 MeV in the case of all the BSk functionals) does 
not seem to have much impact. 
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FIG. 3: Optimal value of nucleon number A as a function of density n at zero temperature in inner 
crust; proton number Z everywhere takes optimal value of 40 for all four forces. 




FIG. 4: Internal energy e per nucleon at zero temperature as a function of density n in inner crust. 
The solid symbols represent the transition densities n^.*^ (see Section IIII CP . 
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FIG. 5: Pressure P per nucleon at zero temperature as a function of density n in inner crust. The 
solid symbols represent the transition densities n^.*^ (see Section IIII C|) . 

B. Continuity with outer crust. 

Our inner-crust code, as used here, is in principle applicable to the outer crust, with the 
background densities nsq vanishing automatically on minimizing the energy per nucleon, 
and it is thus meaningful to compare this code with the code we used for the outer-crust 



calculation of Ref. [7|. In Tabled] we make this comparison at the drip-point density fidrip 
(as determined by the code for the outer crust) with the results for the outer-crust code 
shown in parentheses. 

We see that the inner-crust code (TETFSI) underbinds with respect to the outer-crust 
code (HFB) by around 5 %. This disagreement can be accounted for by the several approx- 
imations made in our TETFSI method, relative to the HFB method adopted in our outer- 
crust calculations [i]], as follows, i) The kinetic energy and spin currents are calculated with 
the semiclassical (T)ETF method, ii) Proton shell corrections are put in perturbatively, and 
neutron shell corrections (shown to be much smaller than proton shell corrections as soon 



as neutron drip sets in [13|, [23(, but obviously not zero, in the outer crust) are neglected 



completely, iii) Rather than allowing arbitrary density variations when minimizing the total 
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TABLE I: Comparison of inner-crust and outer-crust codes at drip point; results for latter code in 
parentheses, e is the internal energy per nucleon, and P the pressure. 



Force 


ridrip (fm 3 ) 


Z 


N 


e (MeV) 


P (MeV fm" 3 ) 


BSkl9 


2.63464 xl(T 4 


40 (38) 


96 (88) 


-1.79426 (-1.87464) 


5.072xl0~ 4 (4.938xl0~ 4 ) 


BSk20 


2.62873xl0~ 4 


40 (38) 


95 (88) 


-1.79451 (-1.87305) 


5.064xl0~ 4 (4.923xl0~ 4 ) 


BSK21 


2.57541 xl(T 4 


40 (38) 


94 (86) 


-1.81718 (-1.90057) 


4.984xl0~ 4 (4.894xl0 -4 ) 


SLy4 


2.45897xl0~ 4 


40 (38) 


93 (82) 


-1.78801 (-1.95898) 


4.744xl0" 4 (4.807xl0~ 4 ) 



energy, the density is parametrized according to Eqs. fl2]) and ([3]). iv) Pairing is neglected 
completely. We have checked that the assumption of sphericity in the inner-crust code has 
a negligible impact in this region of the nuclear chart. 

It will also be seen from Table [I] that there is a slight disagreement in the values of Z 
and N at the drip point. One might speculate that the favoring of Z = 40 over 38 is the 
result of an exaggerated shell effect, but if we drop the proton shell corrections altogether 
then we find slightly higher values of Z, typically 41. However, we have already remarked 
how the inclusion of pairing might well shift the unique value of Z (at T = 0) away from 40, 
and we see from Fig. |2] that a priori it would be difficult to rule out any value of Z between 
36 and 50 at the drip density. The disagreement in the neutron number N is somewhat 
larger, presumably because of our neglect of neutron shell effects, but it is Z that is the 
more astrophysically relevant nucleonic number. 



C. Transition to homogeneous core 



The densities n^.*^f shown in Figs. 0] and |5l and tabulated in Table [TTJ are the densities 



below which homogeneous beta-equilibrated N*M is calculated, for the functional in ques- 
tion, to be unstable to breakup into inhomogeneities. Our values for n^.*^ were calculated 



by the method described in Ref. 33J, in which one defines a free-energy curvature matrix by 



C 



NMe,dyn 



/ 


9fJ, n 






dn„ 


dn p 




dfi p 


d^ P 




dn„ 


drip 


V 









\ 







dn e J 



^ 2Cn n ZC^p ^ 

2C^ n 2Cj p 




4vr 2 e 2 









1 -1 
-1 1 I 



(6) 
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where the fii(= ^-) are the neutron, proton and electron chemical potentials; note that 
n e = n p . The coefficients account for the density-gradient terms in the nuclear density 
functional, which come into play in the presence of inhomogeneities. The third term on the 
right-hand side of this equation gives the Coulomb contribution. Stability of N*M against 
breakup (actually, against density fluctuations of infmitesimally small amplitude) will be 
assured as long as the curvature matrix CNMe,dyn has no negative eigenvalues for all real 
values of k, the wavenumber of density fluctuations. Thus in practice one calculates the 
lowest eigenvalue of CNMe,d y n along the line of /3-equilibrium of N*M in the n n — n p plane 
and determines the density n^.*^ at which it changes sign. Along with n^aaf 5 Table [Til also 
shows the value of the proton fraction Y e and the pressure at the transition point. 

It is instructive to see how our density distributions, as given by Eq. (j2J), approach 
homogeneity as the density increases. In Fig. O we follow the approach to homogeneity by 
showing the neutron and proton density profiles within the WS cell for different values of 
the mean density n. As far as can be seen from this figure, the transition to homogeneous 
matter is very smooth, with no evidence of any discontinuity. However, it is not clear 
in this figure at what precise density homogeneity can be said to set in, but Figs. [7J and 
[HJ complete the picture in this respect. The former shows the variation of the "cluster 
strength" parameters n^ n and function of density: Eq. (121) shows that homogeneity 

corresponds to these parameters being equal to zero. Now in Fig. [7J we see that, for all 
functional, these parameters vanish when the density is very close to n^.*^, calculated 
as described above by a completely different method. A similar conclusion can be drawn 
from Fig. [HJ where we plot a more global measure of the departure from homogeneity, the 
"inhomogeneity factor" 



where V^ e u is the volume of the WS cell and the integration goes over the cell. We plot this 
as a function of density in Fig. [HJ where the transition to homogeneity at a density very 
close to the density n^.*^ is again apparent. 

We stress also that in both Figs. [7J and [HJ the fall to zero of the appropriate measure 
of inhomogeneity is smooth, with no indication of any discontinuity. We cannot exclude 
the possibility that the transition is first order, albeit very weak, but all our results are 
consistent with the transition being of second order or higher. 
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FIG. 6: Profiles of neutron (solid curves) and proton (dashed curves) density distributions in the 
Wigner-Seitz cell for functional BSk21 and different values of the mean density n. Shading denotes 
the region beyond the cell radius. 

Figs. [6] and [7] make it clear that for none of the four functional considered here have 
we found a spherical bubble configuration anywhere in the inner crust. That is, energy 
minimisation always leads to a droplet configuration until homogeneity is reached. This 
result confirms, as far as force SLy4 is concerned, the CLDM calculations of Douchin and 



Haensel [34]. However, we cannot exclude the possibility of very shallow bubbles in a 
very narrow density range, although such configurations would be of limited astrophysical 
interest. Note, moreover, that since our calculations are limited to spherical configurations 
we can say nothing about non-spherical bubbles, such as the very shallow ones found for 
SLy4 in Ref. [16 1. 
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TABLE II: Parameters relating to the crust-core transition. 



Force 


"ffiESf (fm- 3 ) 


-ftrans 


(MeV fm- 3 ) 


BSkl9 


0.0885 


0.0376 


0.428 


BSk20 


0.0854 


0.0356 


0.365 


BSk21 


0.0809 


0.0335 


0.268 


SLy4 


0.0798 


0.0358 


0.361 



IV. DISTRIBUTION OF MASS 



With the EOS determined (for a given functional), the distribution of mass within a 
neutron star (assumed to be non- rotating) is given by the solution to the TOV equations [24 , 



dP(r) Gp{r)M{r) 



dr 



1 + 



P(r) 
c 2 p{r) 



1 + 



4irP(r)r 3 
c 2 M(r) 



1 - 



2 GM{r) 

c 2 r 



and 



Mir) 



47r / p(r')r' 2 dr' 
'o 



Here p(r) is the mass-energy density at the radial coordinate r, given by 

p{r) = n(r) (m + J 



(9) 



(10) 



where M is the nucleon mass and e is the internal energy per nucleon, as plotted in Fig. HJ 
The pressure P{r) appearing in Eq. ([8]) has to be expressed in terms of p(r) through the 
EOS. 

Proceeding as in Section IIIC of Ref. [7], the TOV equations (IE]) and ([9]) are solved for 
the functions p(r) and M.(r) by integrating inwards from the surface (if we had followed the 
usual procedure of integrating outwards from the center our crust results would have been 
contaminated by the uncertainties in the EOS of the core). Then the total baryonic mass 
of the shell of inner radius r and outer radius R, the radius of the star, is 



AM B (r)=4irM r' 2 <S>(r') 1/2 n{r')dr' 

J r 



where we have introduced the metric function 

2GM(r) 



$0] 



1 



c 2 r 



(12) 
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Note that AMe(r), as defined by Eq. (TTTj) . contains the baryonic mass o: 



35 



the entire outer 



36| for SLy4. 



crust, as calculated in Ref. tJ for the three BSk forces, and from Refs. 

We plot AMg(n(r)) as a function of the density n in Fig. Ofor a neutron star of mass 
1.5 Mq and radius 13 km; the fraction of this mass that consists of protons can be read off 
from Fig. [3J given that everywhere we have Z = 40. 

For many purposes it might be more convenient to express AMb as a function of the 
proper depth, given by (see Section 5.6 of Ref. \^\ ) 

which is the only measurable depth in the gravitationally distorted metric. We plot in Fig. ITDl 
n as a function of z, again for a neutron star of mass 1.5 M and radius 13 km, whence 
AMfi(r) can be read off from Fig. [9] as a function of z. 

In Fig. [TT] we show how the total gravitational mass of the crust (inner plus outer) varies 
as a function of the total star mass for a given radius of 9 km. Fig. [12] shows the same 
function for stars of radius 14 km. 



V. CONCLUSIONS 

We have calculated the composition and EOS of the inner crust of neutron stars for the 
three generalized Skyrme-type functionals, BSkl9, BSk20 and BSk21, and for the conven- 
tional Skyrme functional SLy4, using in all cases the TETFSI method at temperature T = 
0. We have also solved the TOV equations to calculate the distribution of mass within the 
crust. 

Qualitatively similar results are obtained for all four forces. In particular, in all cases we 
find Z = 40 for the optimal number of protons per Wigner-Seitz cell throughout the inner 
crust. However, other values of Z lie very close in energy, and if we took pairing into account 
the optimal value of Z might very well be shifted away from 40. Moreover, it is clear that at 
realistic values of the temperature an appreciable range of values of Z will be found. This 
underlines the importance of extending the present calculations to finite temperatures and 
to include pairing. 

The fact that there are no substantial differences in the inner-crust properties for force 
SLy4 and for the three BSk forces despite their having been fitted to different values of 
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the symmetry coefficient J means that this parameter is not of any great relevance in this 
respect. 

We have studied in some detail the transition between the inner crust and the homo- 
geneous core, considering two different measures of the inhomogeneity of our density dis- 
tributions. We find for each of the four functionals that homogeneity is established in our 
calculated distributions at a density very close to the value predicted for the onset in homo- 
geneous N*M of instability against density fluctuations of infinitesimally small amplitude. 

No evidence for bubbles was found in the course of this study of the transition region, 
despite a thorough search. This conclusion does not preclude the existence of non-spherical 
pasta configurations, a possibility that lies beyond the scope of the present paper. Even 
though such phases would have a negligeable impact on the EOS, they might affect transport 
properties. 

The calculations on the inner crust presented here show that our forces BSkl9, BSk20 
and BSk21 make possible a unified and realistic treatment of all regions of neutron stars, as 
in Ref. \M- 




Appendix A: Minimization of Gibbs or Helmholtz functions? 

For simple systems, which in the present context means systems with a single (N, Z) 
configuration, minimizing the Gibbs free energy per nucleon g at constant pressure P is 
completely equivalent to minimizing the Helmholtz free energy per nucleon / at constant 
density n, since in that case the thermodynamic identity 



holds, X denoting any thermodynamical variable. But when two different phases or compo- 
nents, i.e., two different (N, Z) configurations in the present context, coexist in equilibrium 
this identity breaks down, and it is the Gibbs prescription that leads to a correct description 
of the phase transition: there is a discontinuity in the range of densities over which single- 
phase solutions can be found, but the pressure remains constant over this discontinuity, 
which corresponds to the equilibrium coexistence of the two phases. If on the other hand 
one minimizes / at constant density n, discontinuities in the pressure will be found in the 
vicinity of transitions from one (N, Z) configuration to another. An example of this is seen 





(Al) 
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in Fig. [13j where we show the transition from Z = 40 to Z = 20 for functional BSkl4 201 ]. 
with N being optimized in each case. 

Such discontinuities in the pressure are unphysical, and arise in our calculations only 
because our model does not allow the coexistence of two different (N, Z) configurations that 
can occur in reality. But even then, when minimizing / at constant density n, the correct 
equilibrium pressure can be found by making a Maxwell construction, as indicated in Fig. [TBI 
However, on the pressure scale of Fig. [5] these discontinuities will be imperceptible, and the 
Maxwell construction is quite unnecessary: the attendant error will be far smaller than the 
differences between the EOSs of the different functionals seen in Fig. 

In any case, the question of transitions between different values of Z does not arise 
with functionals BSkl9 - 21, since for all these forces Z retains the constant value of 40 
throughout the inner crust. As for changes in N, we recall that this varies continuously in our 
calculations, whence it follows that minimizing / at constant density n leads to absolutely 
no error at all in this respect. 

Appendix B: Pressure formula 

The pressure P at any given point in the neutron-star crust, as given by the EOS and as 
used in the TOV equations, is defined thermodynamically by considering a region of volume 
V that contains the point in question, and is macroscopically sized but small enough for all 
intensive thermodynamic variables to be sensibly constant over it. If F denotes the total 
Helmholtz free energy contained in this region then 

■ < B1 > 

where T is the temperature (here T = 0), N e its number of electrons and N q its number 
of nucleons of type q = n,p for neutrons and protons, respectively. Treating the crust as 
a perfect crystal, this expression remains exact if the region of volume V is taken as the 
appropriate Wigner-Seitz cell, because of the translational symmetry. In the approximation 
used here of spherical WS cells we then have 

P = "4^(m) tm ' (B2) 
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where R is the cell radius. We assume that the Helmholtz free energy in the cell can be 
written as 

r R 

F = 4tt / &rr 2 F{r) , (B3) 
Jo 

where J~(r) is a functional of the nucleon density n q (r) and of the electron density n e (r). 
These densities are related to the total numbers of nucleons and electrons in the cell by 

f R 

N q = 4n / drr 2 n q (r) , (B4a) 



N e = 4n / drr 2 n e {r) . (B4b) 
Jo 

Combining Eqs. (|B2]) and (JB3]) yields 

where 5F/5n q (r) and SF/5n e (r) denote the functional derivatives of F with respect to the 
nucleon and electron densities, respectively. 

Minimizing now the Helmholtz free energy F with respect to arbitrary variations in n q (r) 
and n e (r) leads to the Euler-Lagrange equations 

5F 

5n q (r) 

and 



X q = ^ (B6a) 



5F 

A e = ^ VT , (B6b) 
dn e (r) 

where the X q and A e are Lagrange multipliers introduced to ensure that the nucleon and 
electron numbers given by Eqs. flB4al) and flB4bl) remain fixed; they are identified with the 
corresponding chemical potentials. Using next the identities 

R drr ^^M = ^ nq{R ) , (B7a) 



o 



OR 



and 



R drr 2 ^l = -R 2 n e (R) , (B7b) 
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which follow from the differentiation of Eqs. (1B4aj) and (IB4b[) . respectively, we arrive at 



P = -T(R) + X e n e (R) + \ n z{R) ■ ( B8 ) 

<? 

This pressure formula is a generalization of the expression derived in atomic physics in the 

I I 

framework of the Thomas- Fermi-Dirac model (see, e.g., Ref. [38] and references therein). 

We decompose now the total Helmholtz free energy density in the WS cell into a nuclear 
part, a purely kinetic electron part and a Coulomb part, 

Hr) = ^nuc(r) + J 7 e (r) + Jcoui(r) . (B9) 

Then substituting Eq. (El]) into Eq. flB8]) leads to 

P = -F nuc (R)-F e -F Coul (R) + \ e n e + Y,K n <i( R ) > ( B1 °) 

q 

where we are assuming that n e and T e are position- independent. We now examine in more 
detail the different components of T{f) appearing in Eq. (lB9p . 

In the fourth-order ETF method with Skyrme functionals the nuclear part J r n uc( r ) is a 
local functional of the nucleon densities n q (r) and their derivatives up to just the second 
order, provided the higher-order terms have been integrated as described in Sec. [Ill Note 
that we have not included the proton-proton Coulomb interaction in J r nnc (r). As for the 
electron gas, since it is supposed to be uniform we can write simply 

F e = V7 e {n e ,T) , (Bll) 

where T e is the electron Helmholtz free-energy density, which depends only on the electron 
density n e = N p /V = n p and the temperature T. The Coulomb part of the Helmholtz free 
energy is given by 



R 

2 



^Coul = ^Coul.dir + ■Fcoul.ex = 4tT / .'!/ /' 

Here the direct term is 



J 7 Coul,dir( r ) + J 7 Coul,ex\ r ) 



(B12) 



g 

■Fooni.dirCO = -n c (r)4>(r) , (B13) 

where n c (r) = n p (r)—n e is the net electric-charge density, and 0(r) is the Coulomb potential, 
found on solving Poisson's equation to be given by 

0(r)=47re/ dr' r /2 n c (r')/C(r, r') , (B14) 
Jo 
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in which 

JC(r, r') 

For r = R, Eq. (1B14I) reduces to 



iy j tyt | iy> iy>? | 

2rr' 



(B15) 



<f>(R) 



Aire 



R 



dr r n r (r) = 



(B16) 



the last step being a consequence of global charge neutrality. It then follows from Eq. ( 1B13I) 
that 



? Coul,dir(-R) — 



(B17) 



For the Coulomb-exchange term we have 

T (\ 3e V 3 
>Coul,ex(r) = — - 

4 V 71 



1/3 



xn p (r) i/3 - ^n 4 / 3 



(B18) 



where x is usually equal to 1 but, as explained in Section HI1 is set equal to zero for the BSk 



forces of this paper; for the electrons we have taken the extreme relativistic expression 
Then 
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J 7 Coul(-R) — _ 



3e 2 



7T 



1/3 



xn^ 3 - \nT 



(B19) 



To proceed we have to evaluate the chemical potentials appearing in Eq. ( 1B10I) . The 
Euler-Lagrange equation ( 1B6al) for nucleons can be written explicitly as 



dS/n q {r] 



+ V' 



dV 2 n q {r) 



+ 



xe 



o X 1/3 



S g , P .(B20) 



dn q (r) 

The constant X q can be evaluated at any point r < R, but taking r = R leads to a con- 
siderable simplification of the right-hand side of Eq. (1B20I) . since with our parametrization 
all derivatives of the density vanish at that point. Thus the second and third terms of this 
expression likewise vanish at that point, since each can be expresssed as a sum of terms 
every one of which contains a factor of some derivative of n q (r). Using then Eq. (1B16I) the 
nucleon chemical potential becomes 



9T nuc (R) 
dn q (R) 



— xe 



7T 



1/3 



(B21) 



A further consequence of the vanishing of the derivatives of n q (r) at r = R is that the first 



term here, like the term J-" nuc (i?) appearing in Eq. fIBlOj) . involves only the bulk part of the 
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nuclear free energy density. Next, the Euler-Lagrange equation (1B6bj) for electrons simplifies 
to 



A, 



1/3 



1/3 



(B22) 



because of the uniformity of the electron gas. For the same reason we can write the electron 
pressure (without the Coulomb exchange term) as 



-T e + n, 



(B23) 



dV " ' " dn e 

Also, the Coulomb-potential term e0(r) in Eq. (1B22I) vanishes at r = R, and must be 
negligible for r < R, since otherwise n e and T e would be position-dependent, which would 
be inconsistent with the assumption made and justified in Section [IT] that the electron gas 
is essentially uniform in the inner crust. Then Eq. (1B22j) can be rewritten as 



A e n e = P e + T e + 



e 2 /3 



1/3 



n 



1/3 



(B24) 



2 \tx / 

Substituting now Eqs. flBT9|) . flB2H and flB24l) into Eq. flBTOl) gives us for the total 
pressure 



P — Pnuc + Pe + P( 



Coul,ex 



where 



(B25) 



(B26) 



and 



Coul,ex 



I) 



1/3 



n 



4/3 _ 



x- 



71 



1/3 



nJR) 



4/3 



(B27) 



Given that both terms on the right-hand side of Eq. (IB26I) relate only to bulk matter, 
being independent of any density-gradient terms, it is easy to show from Eq. (1B1I) that P nuc 
represents the purely nuclear pressure of homogeneous nuclear matter with neutron and 
proton densities equal to n n (R) and n p (R), respectively, without any Coulomb term, direct 
or exchange. However, from Eq. (IB18I) it is seen that the last term of Eq. (IB27P is just the 
Coulomb exchange pressure associated with the protons of this homogeneous system, while 
the first term of this equation is likewise the Coulomb exchange pressure of the uniform 
electron gas. 
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This means that the pressure of any crustal layer is the same as that obtained in a homo- 
geneous medium of neutrons, protons and electrons, with the neutron and proton densities 
being those found at the surface of the WS cell, i.e., in the homogeneous background, ng„ 
and n Bp , respectively, while the electron density is to be taken as that of the actual uniform 
electron gas, n e . It is remarkable that the direct Coulomb contribution, calculated exactly, 
vanishes identically, even though n p (R) is not equal to n e . However, this term still mani- 
fests itself indirectly, since it influences the actual values of n n (R) and n p (R) through the 
Euler-Lagrange equations. A similar remark applies also to the inhomogeneities inside the 
cell. 

For the generalized Skyrme force the purely nuclear pressure can be expressed as 

Ti 2 / 5 dC n dC T \ 

Pnuc = xttTq + V C?n 2 m + -C T t n Bt r t + n m — -^n 2 Bt + n B0 —^n Bt r t , (B28) 
3M t t^A 3 dn B0 dn m J 

where n B0 = n Bn + n Bp , while n B1 = n Bn — n Bp , and likewise for r and r 1; with 

r q = ^)^n q {R)^ . (B29) 
The various coefficients are given by 

Cq = |*o + ^sn|o (B30a) 

Cl = ~f° Q + Xo ) ~ 2l t3(1 + X3) ^° (B30b) 
Co = jg*i + (j + x 2 j + —kn B0 + -t s I- + x 5 j nl (B30c) 

CI = ~h Q + + ^ 2 Q + x^j - ^Un B0 Q + x^j + ^t 5 nl Q + z B ) (B30d) 

The pressure P e of the uniform electron gas is calculated as described in Section HI1 using 
expressions given in Section 24 of Cox and Giuli [301 ] . 
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FIG. 7: Variation of the "cluster strength" parameters n\ n and function of density (see 

Eq. ([2]))- The solid symbols represent the transition densities n^.*^ (see Section mrC]) 
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FIG. 8: Variation of the inhomogeneity factor A, given by Eq. (J7]), as a function of density. The 
solid symbols represent the transition densities n^Jf (see Section IIII CP 
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FIG. 9: Variation of baryonic mass of crust (inner plus outer) with density n for neutron star of 
mass 1.5Mq and radius 13 km. 
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FIG. 10: Variation of density n with proper depth z for neutron star of mass 1.5M and radius 13 
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FIG. 11: Variation of gravitational mass of crust (inner plus outer) with total mass of star, radius 
9 km. 
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FIG. 12: As in Fig. \TT\ for stars of radius 14 km. 
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FIG. 13: EOS for functional BSkl4 in the vicinity of the Z = 40 to Z = 20 transition. 
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